Part 2 Fire Selection

2.1 Set Up

2.1.1 Libraries

library(tidyverse)
library(terra)
library(sf)
library(mapview)
library(raster)
library(rgeos)
library(lubridate)
library(ggplot2)
library(exactextractr)
library(patchwoRk)
library(gridExtra)

2.1.2 USDA National Forest Type Group Dataset

Conifer Forest Type Groups: Douglas-Fir, Fir-Spruce-Mountain Hemlock, Lodgepole Pine

# forest type groups and key
conus_forestgroup <- raster('data/forest_type/conus_forestgroup.tif')
forest_codes <- read_csv('data/forest_type/forestgroupcodes.csv')

# set crs
crs = crs(conus_forestgroup)

2.1.3 EPA level-3 Ecoregions

Canadian Rockies, Idaho Batholith, Middle Rockies, Columbian Mountains - Northern Rockies

# level 3 ecoregions
l3eco <- st_read('data/ecoregion/us_eco_l3.shp') %>% 
  st_transform(., crs=crs)

# select northern rocky mountains from level3 ecoregions
eco_select <- l3eco %>% 
  filter(NA_L3NAME %in% c('Canadian Rockies','Columbia Mountains/Northern Rockies','Middle Rockies','Idaho Batholith'))

2.1.4 Mapping

2.1.4.1 Ecoregions

mapview(eco_select)

2.1.4.2 Forest Type Groups

forestgroup_crop <- crop(conus_forestgroup,l3eco)
mapview(forestgroup_crop)

2.2 Define Fire Parameters

2.2.2 Group Adjacent Fires

# function to group adjoining fire polygons to ensure contiguous high-severity patches
group_fires <- function(mtbs_year) {

  # join the polygons with themselves, and remove those that do not join with any besides themselves
  combined<- st_join(mtbs_year, mtbs_year, join=st_is_within_distance, dist = 180, left = TRUE,remove_self = TRUE) %>% 
    drop_na(Event_ID.y)%>% 
    dplyr::select(Event_ID.x,Event_ID.y)
  
  if(nrow(combined)>=1){ # if there are overlaps for this years fires...
    
    # partition data into that that has overlap, and that that does not
    overlap <- mtbs_year %>%
      filter(Event_ID %in% combined$Event_ID.x)
    no_overlap <- mtbs_year %>%
      filter(!(Event_ID %in% combined$Event_ID.x))
    
    print(paste0("there are ",nrow(overlap)," overlapping polygons"))
    
    # join all overlapping features, and buffer to ensure proper grouping
    overlap_union <- st_union(overlap) %>%
      st_buffer(190)
    
    # break apart the joined polygons into their individual groups
    groups <- st_as_sf(st_cast(overlap_union ,to='POLYGON',group_or_split=TRUE)) %>%
      mutate(year = mean(mtbs_year$year),
             Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
      rename(geometry = x)
    
    print(paste0("polygons formed into ",nrow(groups)," groups"))
    
    # join back with original dataset to return to unbuffered geometry
    grouped_overlap <- st_join(overlap,groups,left=TRUE)
    
    # arrange by the new grouping
    joined_overlap_groups <- grouped_overlap %>%
      group_by(Fire_ID) %>%
      tally()%>%
      st_buffer(1) %>%
      dplyr::select(Fire_ID) %>%
      mutate(year = mean(mtbs_year$year))
    
    # add new ID to the freestanding polygons
    no_overlap_groups <- no_overlap %>%
      mutate(Fire_ID = str_c("Fire_",nrow(groups)+c(1:nrow(no_overlap)),"_",year)) %>%
      dplyr::select(Fire_ID,year)
    
    # join the new grouped overlap and the polygons without overlap
    fires_export <- rbind(joined_overlap_groups,no_overlap_groups)
    return(fires_export)
    
    } else { # if there are no overlaps for this year...
      
      print("no overlapping polygons")
      
      fires_export <- mtbs_year %>%
        mutate(Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year)) %>%
        dplyr::select(Fire_ID,year)
      
      return(fires_export)
  }
}
# group adjacent polygons within each fire year
fires_88 <- group_fires(mtbs_select %>%  filter(year == 1988))
## [1] "there are 22 overlapping polygons"
## [1] "polygons formed into 7 groups"
fires_89 <- group_fires(mtbs_select %>%  filter(year == 1989))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_90 <- group_fires(mtbs_select %>%  filter(year == 1990))
## [1] "there are 2 overlapping polygons"
## [1] "polygons formed into 1 groups"
fires_91 <- group_fires(mtbs_select %>%  filter(year == 1991))
## [1] "no overlapping polygons"
# join each fire year, filter by area
mtbs_grouped <- rbind(fires_88,fires_89,fires_90,fires_91)%>%
  mutate(area_ha = as.numeric(st_area(geometry))/10000,
         area_acres = area_ha*2.471)

2.2.3 Select Fires by Ecoregion and Forest Type

# assign ecoregion and proportions of forest type to each fire polygon
fires_join <- st_join(mtbs_grouped,eco_select,join=st_intersects,left=FALSE,largest=TRUE) %>% 
  left_join(., exact_extract(conus_forestgroup,mtbs_grouped, append_cols = TRUE, max_cells_in_memory = 3e+08, 
                             fun = function(value, coverage_fraction) {
                               data.frame(value = value,
                                          frac = coverage_fraction / sum(coverage_fraction)) %>%
                                 group_by(value) %>%
                                 summarize(freq = sum(frac), .groups = 'drop') %>%
                                 pivot_wider(names_from = 'value',
                                             names_prefix = 'freq_',
                                             values_from = 'freq')}) %>%
              mutate(across(starts_with('freq'), replace_na, 0)))
 
# remove unnecessary columns, cleanup names
# filter to ensure fire polygons are at least 25% type of interest
fires <- fires_join %>% 
  dplyr::select("Fire_ID","year","area_ha","area_acres","US_L3NAME","freq_0","freq_200","freq_220","freq_260","freq_280") %>% 
  rename("ecoregion" = "US_L3NAME",
         "freq_df"="freq_200",
         "freq_pp"="freq_220",
         "freq_fs"="freq_260",
         "freq_lpp"="freq_280") %>% 
  mutate(freq_allother = 1-(freq_0 + freq_df+freq_pp+freq_fs+freq_lpp),
         freq_forested = 1- freq_0,
         freq_ideal = freq_df+freq_fs+freq_lpp)%>% 
  mutate(across(starts_with('freq'), round,2))%>% 
  filter(freq_ideal > 0.25)

2.2.4 Select Fires by Burn Severity

# import all mtbs rasters via a list
rastlist <- list.files(path = "data/mtbs", pattern='.tif', all.files=TRUE, full.names=TRUE)
allrasters <- lapply(rastlist, raster)
names(allrasters) <- str_c("y", str_sub(rastlist,22,25))

# create empty dataframe
severity_list <- list()

# loop through mtbs mosasics for 1988-1991
# extract mtbs burn severity raster for all selected fires
# calculate burn severity percentages for each fire
for (i in names(allrasters)){
  mtbs_year <- allrasters[[i]]
  fire_year <- filter(fires, year==str_sub(i,2,5)) 
  raster_extract <- exact_extract(mtbs_year,fire_year, max_cells_in_memory = 3e+09,coverage_area=TRUE)
  names(raster_extract) <- fire_year$Fire_ID 
  
  output_select <- bind_rows(raster_extract, .id = "Fire_ID")%>%
    group_by(Fire_ID , value) %>%
    summarize(total_area = sum(coverage_area)) %>%
    group_by(Fire_ID) %>%
    mutate(proportion = total_area/sum(total_area))%>% 
    dplyr::select("Fire_ID","value","proportion") %>% 
    spread(.,key="value",value = "proportion")
  
  severity_list[[i]] <- output_select
}

# combine extracted raster datasets
severity_df <- do.call(rbind, severity_list) 

# join burn severity % to fires polygons
# fix naming
# filter dataset for 500 acres high severity
fires_severity <- left_join(fires,severity_df,by="Fire_ID")%>% 
  rename(noburn= "1",lowsev = "2", medsev = "3", highsev = "4",regrowth = "5", error = "6") %>% 
  dplyr::select(- "NaN",-"regrowth",-"error") %>% 
  mutate(highsev_acres = area_acres*highsev)%>% 
  filter(highsev_acres > 500)

2.2.5 Clean Up Dataset

# get the most common forest type within each polygon
fires_select <- fires_severity %>%
  left_join(.,exact_extract(conus_forestgroup,fires_severity, 'mode', append_cols = TRUE, max_cells_in_memory = 3e+08)) 

fires_select$mode <- as.factor(fires_select$mode)

fires_select <- fires_select %>% 
    mutate(fire_foresttype = case_when(mode==200 ~ "Douglas-Fir",
                                       mode==220 ~ "Ponderosa",
                                       mode==260 ~ "Fir-Spruce",
                                       mode==280 ~ "Lodegepole Pine",
                                       TRUE ~ "Other"),
           Fire_ID = str_c("Fire_",c(1:nrow(.)),"_",year))
# join the grouped fires back to original mtbs boundaries
fires_mtbs <- st_join(mtbs_select,fires_select,left=FALSE,largest=TRUE) %>% 
  filter(year.x==year.y)%>% 
  dplyr::select("Event_ID","Incid_Name","Fire_ID","Ig_Date","year.y","state","BurnBndAc","ecoregion") %>% 
  rename(year= year.y)

2.3 Final Fire Dataset

2.3.1 Selected fires by year

# plot
mapview(fires_select, zcol = "year")

2.3.2 Selected fires by majority forest type

# plot
mapview(fires_select, zcol = "fire_foresttype")

2.4 Export Data

2.4.1 Final Cleanup for Export

# reformat and project
fires_export <- fires_select %>% 
  mutate(year = as.integer(year)) %>% 
  st_transform(., crs="EPSG:4326")

mtbs_export <- fires_mtbs %>% 
  mutate(year = as.integer(year)) %>% 
  st_transform(., crs="EPSG:4326")

2.4.2 Export

# st_write(fires_export, "data/fire_boundaries/", "fires_export.shp", driver = 'ESRI Shapefile')
# st_write(mtbs_export, "data/fire_boundaries/", "mbts_export.shp", driver = 'ESRI Shapefile')